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Differences in masses inferred from dynamics, such as velocity dispersions or X-rays, and tliose 
inferred from lensing are a generic prediction of modified gravity theories. Viable models however 
must include some non-linear mechanism to restore General Relativity (GR) in dense environments, 
which is necessary to pass Solar System constraints on precisely these deviations. In this paper, we 
study the dynamics within virialized structures in the context of two modified gravity models, f{R) 
gravity and DGP. The non-linear mechanisms to restore GR, which f{R) and DGP implement in 
very different ways, have a strong impact on the dynamics in bound objects; they leave distinctive 
signatures in the dynamical mass-lensing mass relation as a function of mass and radius. We present 
measurements from N-body simulations of f{R) and DGP, as well as semi-analytical models which 
match the simulation results to surprising accuracy in both cases. The semi-analytical models are 
useful for making the connection to observations. Our results confirm that the environment- and 
scale-dependence of the modified gravity effects have to be taken into account when confronting 
gravity theories with observations of dynamics in galaxies and clusters. 
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I. INTRODUCTION 

Gravity, as described by General Relativity (GR), is re- 
markably weakly constrained in the present day on scales 
larger than a few AU. Though measurements from bi- 
nary pulsar timing to the cosmic microwave background 
(CMB) and big bang nucleosynthesis are all consistent 
with GR, there is still room for order unity deviations 
in the cosmos today, on scales of kpc and larger. Thus, 
testing gravity on cosmological scales is an interesting 
frontier, and the focus of much current research [l|, [2, 




' Any gravity theory that attempts to be complete has 
to satisfy stringent Solar System constraints, and has to 

] locally match the predictions of GR to within one part 
in 10^ there. Only a few consistent models which modify 
gravity appreciably on large scales, but restore GR lo- 
cally are known. Two of them will be the subject of this 

■ study: f{R) gravity [IMl, and the DGP model 

', Within certain bounds placed by the CMB and expan- 
sion history measurements in addition to Solar System 

\ tests, both theories can be made to satisfy all current 
constraints on gravity (including the observation of an 
accelerating expansion). In both models there exists a 
non-linear mechanism to restore GR in high-density en- 
vironments: the chameleon effect for f{R), and the Vain- 
shtein mechanism for DGP. Furthermore, all currently 
known consistent modifications of gravity on large scales 
include some variant of either of these mechanisms. In 
order to be able to constrain these models with cosmolog- 
ical data, it is crucial to correctly include the non-linear 
mechanisms. Recently, N-body simulations of /(i?) [13 
and DGP [l8l - l20j have been done which self-consistently 
solve the non-linear field equations together with the 
growth of structure (see also [2l| for the first study of 
the DGP case, using a different approach). In principle. 



it has become possible with these simulations to unlock 
the wealth of observations available on non-linear scales 
to probe gravity, albeit in a necessarily mo del- dependent 
way. 

It is well known that the additional degrees of freedom 
present in modified gravity theories generically affect the 
dynamical potential, which governs the propagation of 
non-relativistic bodies, differently than the lensing po- 
tential, which governs the propagation of massless parti- 
cles such as light (e.g., ^^). Thus, comparing dynamical 
with lensing mass estimates is an interesting and quite 
generic probe of modifications to gravity. In this paper, 
we study the signatures of f{R) and DGP in dynami- 
cal observables such as velocity dispersions, compared to 
lensing which measures essentially the "true" mass (i.e. 
the integral over the rest-frame density) in both models. 

Constraints on the difference between dynamical and 
lensing potential are often phrased in terms of the post- 
Newtonian parameter 7ppN [Eq. ([5]) below] , in analogy to 
Solar System tests. In general, however, the departures 
from GR cannot be encapsulated by a single parame- 
ter but are functions of scale, time and the local envi- 
ronment. In particular, this is the case for both f{R) 
and DGP. Hence, we introduce a more generally applica- 
ble quantity ^ [Eq. ([3])] which is defined directly via the 
modified forces, and is well suited for predictions in the 
context of f{R) and DGP as well as for constraints from 
observations. 

Velocities of extragalactic objects are measured 
through their redshifts z, which receive a contribution 
|Az| = "Uii/c from the line-of-sight velocity In the 
cosmological context, there are two regimes where the 
dynamics of matter can be understood fairly easily: on 
very large scales, linear perturbation theory in the mat- 
ter density is valid, simplifying the theoretial predictions. 
Large-scale velocity fields can be measured through the 
rcdshift distortion of the power spectrum, which thus 
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offers a probe of the dynamical potential [23|, |2J]. On 
small scales, most of the observable matter lies in grav- 
itationally bound dark matter halos. In this regime, for 
relaxed systems, the velocity distribution of coUisionless 
objects such as dark matter, galaxies or stars is related 
to the dynamical potential by the virial theorem. For 
coUisional particles such as diffuse gas, this relation is 
given by hydrostatic equilibrium. The virial or thermal 
velocities can be observed as velocity dispersion of stars 
in galaxies, galaxies in clusters, or as X-ray or Sunyaev- 
Zeldovich signal from diffuse gas in clusters. Also, the 
redshift-space matter pow er sp ectrum on small scales is 
a probe of virial velocities j23l [25|. 

This paper is concerned with the latter regime, and 
our goal is to study the dynamics of matter in halos. 
Since these are highly non-linear systems, rigorous re- 
sults can only be obtained via N-body simulations. We 
therefore present measurements from the modified grav- 
ity simulations of f{R) and DGP [13, [11 HI. However, 
for many practical purposes including comparison with 
observations, it is necessary to go beyond the simulation 
results which have limited resolution and cover only a 
few points in the parameter space of the models. Thus, a 
sufficiently accurate semi-analytic model of the dynam- 
ics in modified gravity is desirable to bridge the gap with 
observations. 

Fortunately, we can make some justified assumptions 
which simplify the problem greatly: first, since we are 
concerned with sub-horizon scales, we employ the quasi- 
static approximation, neglecting time derivatives and as- 
suming the halos are in steady state. Further, we as- 
sume spherically symmetric halos. While certainly not 
realistic, deviations from spherical symmetry are not ex- 
pected to affect the results qualitatively. Throughout, we 
will assume a Navarro-Frenk- White (NFW) [23 profile, 
although all derivations can easily be generalized to dif- 
ferent profiles. The problem is then reduced to finding 
the solution of the field equations for a spherically sym- 
metric mass, and calculating the modified gravitational 
force. The accuracy of this simplified model can then be 
benchmarked with the simulation results. 

The paper is structured as follows. In § Ull we in- 
troduce our main observable, the modified gravitational 
force strength, and present the theoretical expectations 
and semi-analytic models for f{R) and DGP. §|lTT] con- 
tains the simulation results and comparisons with the 
theoretical models. We then discuss the application to 
observations in § IIVI We conclude in § |Vl 



II. THEORETICAL EXPECTATIONS 

In this section, we derive theoretical expectations for 
the modified gravitational forces and virial quantities 
measured in the simulations in § IIIII and connected to 
observations in § IIVI Gravitational forces arc given by 
the gradient of the dynamical potential ^I^, defined via 



the perturbed FRW metric in Newtonian gauge: 

ds^ = -(1 + 21')(ii2 + a^{t){l + 2$)dx2. (1) 

As a reference point, we consider General Relativity (GR) 
in the Newtonian limit, where the dynamical potential 
satisfies the Poisson equation 



V^^- = V^I-AT = 4:TtG5p, 



(2) 



where 6p = p/p is the total matter overdensity. Assum- 
ing spherical symmetry, which we will throughout, we 
can define a parameter g : 



lir) 



d-^/dr 
d^N/dr' 



(3) 



which quantifies the strength of the gravitational force 
in modified gravity relative to that which would be mea- 
sured in GR given the same density field, g = 1 corre- 
sponds to unmodified forces. Here we have suppressed 
the dependence of g on the scale factor a; unless other- 
wise stated; wc will always assume a = 1. 

In the models we consider, the lensing potential satis- 
fies^ 



N- 



(4) 



Hence, g can be probed for example by comparing dy- 
namical to lensing mass estimates of a given object. Such 
cornparisons in the Solar System or for distant galax- 
ies [28| are often phrased in terms of the post-Newtonian 
parameter 7ppn: 



7PPN 



) 



1^=°2 



1. 



(5) 



The last equality relating 7ppN to our g parameter (where 
we have used Eq. Q) is only valid when the force mod- 
ifications are scale-independent, such as in Brans-Dickc 
(BD) type scalar-tensor theories. Note that the PPN pa- 
rameter is formally defined via the potentials, while our 
^ parameter is derived in terms of forces. Only forces, 
or more generally derivatives of the potentials ^f, are 
observable, and a specific solution of the potentials (e.g., 
the Schwarzschild metric) is used to infer 7ppn- However, 
in the models we consider ^ is generally scale-dependent, 
i.e. the scalar degrees of freedom do not follow the same 
scaling with distance as the GR potentials. Hence, it is 
advantageous to define a parameter based directly on the 
forces, rather than 7ppN which is not immediately linked 
to observables. 

In many practical cases, one is interested in a weighted 
average of ^ over an object or region of space. 



J r'^wir) ^(r) dr 
J 7'^w{r) dr ' 



(6) 



^ In f{R), there are corrections of order f^, and < \ffto\ < 
10~* for the models we consider; this is negligible compared to 
the 0(0.1) effects we will discuss. 
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where w is a weighting function depending on the pre- 
cise observable considered (we will turn to this in ^ lIVp . 
The key point is that given a prediction for ^(r) we can 
estimate any such weighted average (as long as spherical 
symmetry is a sufficiently good approximation). In the 
next section, we will introduce one such averaged force 
modification which is relevant for comparison with sim- 
ulations. Wc will then review the Newtonian potential 
and scaling relations for a dark matter halo with NFW 
profile, before studying the same case for f{R) and DGP 
gravity. 



A. Virial theorem and velocity dispersion 

For the comparison with our dark matter-only simu- 
lations, it is useful to consider a coUisionless system in 
virial equilibrium. In that case, the virial theorem states 
that 



W 
W 

T 



-2T, where 



- J d^x p(x) X • V'I'(x), 



(i3xp(x)cr^3£,(x). 



(7) 

(8) 

(9) 



denote the trace of the potential energy tensor and po- 
tential energy, respectively. Here = Scr^ is the 
three-dimensional velocity dispersion (see Eq. ([57)) for 
our practical definition in terms of dark matter parti- 
cles). Since the virial theorem is derived from the coUi- 
sionless Boltzmann equation, and is thus a consequence 
of energy-momentum conservation, it is unchanged in any 
metric theory of gravity, and hence also in the models we 
consider. The modification enters through the modified 
relation between the potential and the matter distri- 
bution. 

Note that in the cosmological context, we are not deal- 
ing with strictly isolated systems, so that Eq. ([7]) does not 
hold precisely. Nevertheless, the validity oi W = aT for 
simulated dark matter halos has been shown to hold to 
high accuracy [29l - [3l| . Here, the constant a depends on 
the mass and radius definition chosen for the halos. 

In the spherically symmetric case, we can use the defi- 
nition of ^ [Eq. ([3])] to relate the potential energy tensor 
and kinetic energy in modified gravity to the Newtonian 
values Wn, T^- 



(10) 



where J^.^. is given by Eq. (jS]) with a weighting function 



Wvirir) = p{r) r- 



dr 
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The gradient of the Newtonian potential appearing here 
is uniquely determined by the density p(r), assuming that 
external tidal fields are negligible. 



B. Newtonian potential of a halo 

Let us consider the GR case first. We can integrate 
Eq. © to obtain: 



d^N _ GdM{<r) 



dr 
SM{< r) 



47r 



dr' r'2 Spir'). 



(12) 
(13) 



Note that SM is defined in terms of the enclosed over- 
density dp. Imposing the condition 5'iv('' oo) = 0, we 
can integrate again and obtain: 



*Ar(r) 



dr 



,G5M{< r') 



(14) 



Let us now consider an NFW halo with mass Ma, de- 
fined as the mass contained within a radius i?A so that 
the average density within i?A is p A (note that A here 
is arbitrary and does not have to correspond to a certain 
"virial" overdensity) . The NFW profile has been shown 
to be a good match even to the halos in modified gravity 
simulations [lO,!!^- We define the corresponding concen- 
tration as ca = Ra/ts- We will consider an untruncated 
profile here; while this overestimates the exterior mass 
somewhat, it is closer to the profiles measured in simu- 
lations than the other simple choice, a truncated profile. 
Then, the density profile is given by 



p(r) 
/nfw(2/) 



4p,s Infw {r/rs), 
1 



(15) 
(16) 



where ps = pifs) is chosen so that the mass within i?A 
is AfA, and we have 



5M{< r) = A/a 

P{y) = -T- 



FjcAr/RA) 
ln(l + 



1 - 



F{ca) {r/R^f 



F{cAr/RA) A 



y 



(17) 
(18) 



The correction in square brackets in Eq. (|17p is usually 
neglected since it is smaller than A~^, and we will do 
so here as well in order to simplify the analytical expres- 
sions. From this, we get 



d^ 



N 



r2 



dr Ra 7-2 F{ca) ' 
where the potential scales with ^I'a defined by 

GMa 



^'a 



Ra 



(19) 



(20) 



We can integrate to obtain the potential for an isolated 
NFW halo: 



*Ar(r) 

E{x,c) 



-^aE 

(1 + 



c) ln(l -I- cx) 



(1 + c)x ln(l + c) — cx 



(21) 
(22) 
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£■(0, c) 5- 12 (for c 4-30) gives the central depth of 
the potential well for an isolated NFW halo with concen- 
tration c in units of ^a- Note that in reality, the depth 
of the potential well will depend on the large scale envi- 
ronment, so that Eq. (|2ip will only give a rough scaling. 
^'a in turn is given by: 



*A = {GM^Hof^ ( |anA 



1 



1/3 



Ma 



2/3 



6.26 X 1022 A'fo//i 
A/a 



1.79 X 10" 



1015 Afo//i 



2 A 

2/3 



1/3 



(23) 

(24) 
(25) 



where for the last equahty we have assumed = 0.25 

2/3 

and A = 200. Note the scaling of ^'a with , because 
Ma and Ra are linked through the fixed overdensity A. 
Throughout, unless otherwise stated we use the concen- 
tration relation of 1331: 



c(M, z) 



9 



1 



M 



-0.13 



(26) 



Here, M* sa 3.2 x 10^2 Mq/H for our fiducial ACDM cos- 
mology. Recently, more accurate expressions for the con- 
centration have been found (s^ . [3a |. However, our re- 
sults are not very sensitive to the concentration, hence 
we deem Eq. ((25)) sufficient. At the very highest masses 
Ma > 10^5 Mq/H however, Eq. pB]) underpredicts the 
concentration significantly (e.g., [HjIsBl)- As a simple 
remedy, we take c = max{4, c(M)} in place of c{M) from 
Eq. (Ell). 

Finally, the weighted g{r) quantifying the modification 
to the potential and kinetic energy [Eq. ([T0| ] can be writ- 
ten as 

_ _ Jq dx X F{cAx)f]^Fw{cAx)g{x Ra) ^^y-j 
Jq dx X F{c ax) f NFW (cAx) 



where 



A- 



C. f(R) 

f{R) gravity (see [lH| for a review) is a modified action 
theory where the Einstein- Hilbert Lagrangian R/IQttG is 
replaced with [R+ f{R)]/167TG. Throughout this section 
R denotes the Ricci scalar. f{R) models correspond to 
scalar-tensor theories, where the scalar degree of freedom 
is given by = df/dR and mediates the relation be- 
tween density and space-time curvature. In order for the 
theory to be stable under perturbations, it is necessary 
that /ii < [13. 

In the smooth background, the scalar field assumes a 
value of f]i = fR{R); where R oc is the scalar curva- 
ture of the background. In this paper, we only consider 
models with < 10"'^, and will thus drop higher order 



terms in the fa field which simplifies the expressions. In 
the quasi-static regime, the /r field and the dynamical 
potential are then determined from the density field by 
the following coupled equations: 



^^SfR 



1 



[SRUr) - SwGSp] , 



V^* = ^-^^Sp^lsRifR). 



(28) 
(29) 



Here, S stands for perturbations from the background 
value: SJr, = /r — Jr and 6R = R — R. R and SR, are 
non-linear functions of the field Jr, hence Eqs. (|28 |) - (p9)) 
are difficult to solve in general. However, there are two 
limiting cases which can be solved easily. 

First, consider the case where Jr is much larger than 
typical potential wells in the universe. In that case, SJr 
sourced by the r.h.s. of Eq. ([25)1 is always much less than 
/r, and we can linearize the 6R term: 



1 



fRRiR) 



SfR, 



(30) 



where Jrr = d^f/dR"^. Eq. ([28| then becomes an equa- 



tion for a massive scalar field with m 



We call the inverse mass Ac of the field in the background 
the Compton wavelength. In this limit, 6R ^ Ron scales 
smaller than Ac. Eq. (gH]) then tells us that * = 4/3\l>jv, 
i.e. gravitational forces are increased by 4/3 within the 
range of the Jr field given by Ac. 

In the opposite limit, both terms on the r.h.s. of 
Eq. (|28p are much larger than the l.h.s. ~ fR./r'^ on 
the scales of interest. Since the field perturbation is lim- 
ited in magnitude to be less than 5fR has to adjust 
itself so that the two terms on the r.h.s. cancel to a high 
degree, in other words 



5R{fR) ~ 8ttG6p. 



(31) 



Hence, the GR expression is restored, and Eq. yields 
^' = AT accordingly. This is called the chameleon regime 

H. 

In order to determine the transition between these two 
regimes, we consider the solution for a spherically sym- 
metric mass. Formally, we can write the solution for S/r 
as: 



SfR{r) 



2 GSM,s{< r) 



6M,ff{r) = 4tt / dr' r'^5p,s{r'), 
Jo 

6R{r) 



Spcs{r) = Sp{r) - 



SttG 



(32) 
(33) 
(34) 



With these definitions, the modified dynamical potential 
satisfies 



(35) 



Eqs. (IMll-dMl) state that SM^g < SM. If the pertur- 
bation SfR is small for all r (which in general is only 
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true far away from the body), we can neglect the SR 
term in Eq. ([M]). Then, SM^g = SM and we have 
I^/-r('')I " S/Sj^E'jvl- However, the maximal value Sfu 
can achieve is in which case the fa field becomes 0. 
Thus, we arrive at the following condition: 

|/i?o| < l^N. (36) 

If the value of ^'at for the body is larger than this, the 
field must enter the chameleon regime. Then, Spcs is non- 
zero only outside of the radius where Eq. p6p is met. 
(5Meff is thus given by the mass outside of this radius 
which can be thought of as forming a thin shell. For this 
reason, Eq. is also called the thin-shell criterion. 

Since cosmological potentials range from 10~^ — 10~^, 
we expect that the chameleon mechanism will operate 
for background field values < 10^^. 

This general picture holds for any viable functional 
form of f{R)- However, in order to evaluate the effect 
on the dynamics quantitatively and to compare with the 
N-body simulations, we have to adopt a specific model. 
The functional form used in the simulations [13, IH, 
is the one of [s^l with n = 1, i.e. 

parametrized by the two constants A and i?c- If the 
present-day background curvature is much greater 
than Rc, which will be the case for the f{R) models con- 
sidered here, we can expand Eq. (j37p to first order in 
Rc/R, and define a new parameter fnQ — f{Ro) so that 

/(i?) = -2A-/;^o^. (38) 

The first term supplies an effective cosmological constant 
yielding accelerated expansion of the background. The 
second term, controlled by fuo <^ 1 determines the de- 
partures from GR, and yields corrections to the back- 
ground expansion of order fnQ. Since we will consider 
models with \fRo\ < 10^^, the background expansion 
is essentially indistuingishable from ACDM. Taking the 
derivative of Eq. ([55)1 . we obtain the relation between the 
scalar field and the local curvature at the present day: 

fR = /™||. (39) 

Furthermore, the Compton wavelength Ac oc \/ Jrr (x 

Simulations were performed for a range of background 
field values |/i^o| = 10"^, 10-^ 10""*. From our discus- 
sion above, we expect that the chameleon mechanism will 
operate in the intermediate and small field cases, while 
it will be essentially absent for the large field (10"'*). In 
addition to the f{R) simulations, ordinary ACDM simu- 
lations were performed using the same expansion history 
and initial conditions. The cosmological parameters used 
in the simulations are summarized in Tab. ID 




r/R 



FIG. 1: Spcs [Eq. (|34p ] divided by the mean matter density 
p, determined from the numerical solution of the f{R) field 
equation for an NFW halo with M300 = 2 x 10^"^ Mq /h for 
different values of fm. Also shown is the matter density ^p/p 
of the halo itself (dotted line almost matching the 10""* field 
curve). The dash-dotted line shows a density profile which 
matches that measured in simulations, including an additional 
external overdensity [Eq. (|4T]) ] . 

Given a density field such as that for an isolated NFW 
profile, one can solve Eq. p8|) numerically. We have 
done so for the spherically symmetric case using a one- 
dimensional relaxation algorithm (in fact we solve for u 
defined by fR = exp(u) to avoid overshooting to /j^ > 
[l7|). While only an approximation of the physical re- 
ality, the spherically symmetric case allows for a much 
higher resolution (at much smaller computing time) than 
achievable in the full 3D cosmological simulations. The 
boundary conditions are given by dfR/dr = at r = 0, 
and SJr = at the outer edge of the grid, chosen here 
to be Tiiiax = 50Mpc//i. We use 4096 equally spaced 
grid points in r. Once SJr is known, Eqs. can 
be evaluated using SR{SfR), and the modified forces are 
given by 

, , 1 Moff (< r) , , 

Fig. [1] shows the "effective density" (5pcff which sources 
the perturbation 8jR to the field, for a halo of mass 
2 X 10^^Mq//i and different values of ]rq. For large 
values of |/i?o| 2 x 10^^, the thin shell condition is 
never met, so that (^poff = everywhere (except at very 
large r where the field decays due to its finite Ac)- For 
smaller field values, we can see that a "thin shell" devel- 
ops. For |/i?o| = 10~^ it is quite broad, while it narrows 
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FIG. 2: jjf(R) as a function of the scaled radius r/Rsoo for 
|/flo| = 10~^ and different halo masses, from the numeri- 
cal spherically symmetric solution. The low-mass halo is un- 
screeened, showing the 4/3 force enhancement throughout, 
while higher mass halos are partially screened. The arrows 
for the two more massive halos indicate at which r the condi- 
tion Eq. ((36)) is first met. For the 2 x 10^" M^/h halo, we also 
show gf(R) including an external density field (dash-dotted 
line; see Fig. [T] and Eq. (STJ). 



FIG. 3: Averaged modified gravitational force jjvii.f(R) for 
NFW halos as a function of halo mass for different values of 
\fRo\- The points and thick lines show the numerical results 
from the ID relaxation code. The thin black lines show the 
predictions of a simplified model Eq. (|42|) . The behavior with 
mass is similar in the different models, with the transition 
mass between unscreened and screened regimes shifting as 
expected by simple estimates. 



considerably for a small field of |/flo| = 10~^. Note that 
the transition to peff = within the shell is very sharp, 
owing to the much smaller Compton wavelength within 
the body (recall Ac oc i?-^/^). 

Fig. [5] shows gf{R){r) for \fBa\ = 10^^ and differ- 
ent halo masses. The IQ^^ M^/h halo is unscrcccncd, 
showing the 4/3 force enhancement throughout. The 
2 X 10^'^ Mq/K halo is partially screened, while the 
10^^ Mq/H halo is screened to a large extent within i?3oo- 
For the latter two cases, we also indicate the screening 
radius Tgcr where, going from the outside in, the thin 
shell condition Eq. is first met. This radius serves 
as an indication of whether a given mass is screened, and 
roughly to what extent. As Eq. (pS]) shows, the screening 
radius depends on the depth of the potential well, which 
is influenced by the large-scale environment. To inves- 
tigate this effect, we have added an external large-scale 
density field to the NFW profile, roughly matched to the 
halo profiles in our simulations [32| at large radii: 

where 5/Onfw is the halo overdensity given by the NFW 
profile, and the external density is smoothly cut off at 



^ 40Mpc/ft.. (5/9h-fcxt is shown as dash-dotted line in 
Fig.[T] The resulting ^^(fl-) including the external density 
field is shown in Fig. [2] for the intermediate mass halo. 
Clearly, the field is screened at somewhat larger radii in 
this case, and sf{R) is smaller than that predicted for the 
NFW profile alone by about 0.04 in the transition region. 
Since halos can reside in a variety of environments, we 
expect significant scatter in the strength of the modified 
forces within halos in the f{R) case, halos in overdense 
regions being screeend more strongly than those in av- 
erage or underdcnsc regions. Further, we expect that 
the environment-dependence will be more significant for 
lower mass halos than for massive halos (> 10^'^ Mq/H), 
since the former can be affected by a massive halo nearby, 
while the latter usually dominate their environment. We 
also investigated the effect of varying the halo concentra- 
tion by ±20%; the impact on ^/(i?) is small in comparison 
with the effects of the large scale environment however. 

Finally, using the results for ^f(B,}(j') we can evaluate 
Eq. (P7)) . Fig. [3] shows ^vir,/(i?,) as a fimction of mass 
for different values of the background field fno. The 
thick lines and points show the numerical results from 
the relaxation code. For the strongest field, only the 
most massive halos (more massive than found in our lim- 
ited volume simulations) arc chameleon-screened. For 
the weakest field, all halos above AI ^ 10^^ M^/h are 
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expected to be screened, while for the intermediate field 
the transition scale is around lO^"^'^ Mq/H, relevant for 
galaxy clusters. Wc will compare the predictions for both 
£ (r) and J . with simulation results in § IIIII 

As a simple analytic model for the numerical results, 
we make the assumption that all mass of the halo outside 
of Tscr contributes to SMcs- This results in the following 
simple prescription: 



1 M{< r) 
3~ 



M{< r^cr) 



1 



1 



M{< r) 
F{cAr/RA) 



(42) 



We then form the same average via Eq. pT]) . As shown in 
Fig. 131 (thin black lines), this approximation predicts the 
onset of the chameleon screening quite well, though the 
predicted transition between unscreened and screened 
regimes is somewhat too sharp. Nevertheless, this simple 
model can be useful in interpolating the numerical results 
for different values of /aq- 



D. DGP 

In the DGP brancworld scenario [l^ , matter and radi- 
ation live on a four-dimensional branc in five-dimensional 
Minkowski space. The action is constructed so that on 
scales larger than the crossover scale rc, gravity is five- 
dimensional, while it becomes four-dimensional on scales 
smaller than rc- This model admits a homogeneous cos- 
mological solution on the brane which obeys a modifed 
Friedmann equation [40l |: 



H^±—=8ttG[p + pbe]- 
rr 



(43) 



TABLE I: Parameters of the simulated f{R) and DGP cos- 
mologies. For each model, GR simulations with identical ex- 
pansion history and initial conditions were also performed. 







sDGP 


nDGP-1 


nDGP-2 




0.24 


0.258 


0.259 


0.259 


(eff.) 


0.76 





0.741 


0.741 


lg|/fio| 


-4, -5, -6 








rc [Mpc] 




6118 


500 


3000 


/3(a = 1) 




-1.15 


1.21 


2.25 


Ho [km/s/Mpc] 


73.0 


66.0 


71.6 


71.6 


100 ^It h' 


2.23 


2.37 


2.26 


2.26 


Us 


0.958 


0.998 


0.959 


0.959 


lif yla(0.05 Mpc^^) 


2.24 


2.02 


2.11 


2.11 


cr8(ACDM)'' 


0.796 


0.657 


0.789 


0.789 



"Linear power spectrum normalization today of a ACDM model 
with the same primordial normalization. 



The sign on the l.h.s. is determined by the choice of 
embedding of the brane. The negative sign is called the 
self- accelerating branch, since it allows for accelerated 
expansion even in the absence of a cosmological con- 
stant. The positive sign is called the normal branch, 
which does not exhibit self-acceleration. Here, we con- 
sider models of both branches (see [3 US, Ej ) : a self- 
accelerating model without a A term (pde = 0), sDGP, 
where rc ^ 6000 Mpc is adjusted to best match CMB and 
expansion history constraints [i^l (note that this model 
is in ~ 4 — 5(7 conflict with current data) ; and normal- 
branch models with a dark energy component pde ad- 
justed so that the expansion history is exactly ACDM 
20l |. In that case, rc is a free parameter, and we chose 
values of 500 Mpc {nDGP-1) and 3000 Mpc {nDGP-2). 
The remaining cosmological parameters are summarized 
in Tab. [D For both sDGP and uDGP models, wc have 
also performed ordinary GR simulations employing the 
same expansion history and initial conditions as for the 
DGP simulations. 

On sub-horizon scales, and scales smaller than the 
crossover scale r^ DGP braneworld models can be ac- 
curately described as scalar-tensor theory [i^l , where the 
brane-bending mode ip mediates an additional attractive 
(normal branch) or repulsive (self-accelerating branch) 
force. Gravitational forces in DGP arc governed by: 



N ■ 



(44) 



The Lp field is sourccd by matter overdensitics simi- 
larly to the usual GR potentials, but has quadratic self- 
interactions which suppress the field once density con- 
trasts become non-linear. The full equation for the Lp 
field is (assuming a = 1; see e.g. Q): 

SttG 



VV+ ^[(VV)^ - (V,V,^)(V'V»] 



5p. (45) 



Here /3 is determined by the expansion rate H{a) via 

H{a) ^ 



/3(a) ^l±2H{a)rc 1 



3i?2(a) 



(46) 



where the positive (negative) sign are valid for the normal 
(self-accelerating) branch. The present-day values for /3 
are given in Tab. ID 

While analytical solutions to Eq. psj) do not exist in 
the general case, the case of a spherically symmetric mass 
is solvable in terms of closed expressions [3, EH . In par- 
ticular, one obtains the following equation for the gradi- 
ent of Lp 4l|: 

dip _ G5M{< r) 4 
dr r^ 



W \r.{r) 



9(0 = e{Vi+F^-i) 



(47) 
(48) 



r*(r) in Eq. (|T7| is the r-dependent Vainshtein radius 
defined as: 

1/3 



r*{r) = 



r i6GSM{< r)r 



(49) 
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Note that r/r* is a function of the average overdensity 
5p{< r) within r. Specifically, scaling to a halo with mass 
Ma and radius i?A determined by a fixed overdensity A 
and neglecting the small difference between M and 5M , 
we have: 

where x = r/R^ and the quantity e is determined by the 
background cosmology: 

£ -^iHorc)^^rna-^ (51) 

At a = 1, = 0.32 for sDGP, and 0.002/0.023 for nDGP- 
l/nDGP-2, respectively. Using Eq. ©, (gT]), we 

then have 

.DGPW-i + ^.(;r^). (52) 

On large scales where 6p{< r) <^'p, r is much larger than 
so that g{^) — > 1/2 and dip/dr becomes simply pro- 
portional to d'^M/dr. In this limit, ^dgp = <9DGP.iin = 
1 + 1/ (3/3). This is the same expression one would obtain 
by simply neglecting the non- linear terms in Eq. (|45p. 
On small scales where r ^ r*, modified forces are sup- 
pressed by (e^)^^/^, where 5 = 5p{< r)/'p is the average 
overdensity within r. 

Note that the specific tensorial structure of the non- 
linearities in Eq. (|45p is crucial to recover the linearized 
expresion ^DGP,iin- It is possible to simplify Eq. (|45p by 
neglecting the tensorial structure, resulting in a Poisson 
equation for Lp with a source term given by a non-linear 
function of 5p [2l| . However, this simplification quali- 
tatively changes the large-distance behavior of the Vain- 
shtein mechanism [2l1. l4l|. Thus, it will turn out to be 
crucial that the simulations solve the full Eq. (|^5|) for our 
comparison with the theoretical predictions from Eq. (|T7)) 
and Eq. 

Note that in the Vainshtein limit, 

(p{r<^^r^)^C ^ -(x-^Nir) — , (53) 

r* r 

where C is a constant of order unity [4l|. Hence, the 
ip field itself is suppressed less than the modified forces 
by a factor of (r/r*)^/^. However, only the forces are 
observable. This shows that in theories with non-linear 
interactions such as DGP, quantifying departures from 
GR in terms of forces is more appropriate than the pa- 
rameter 7PPN defined in terms of the potentials [Eq. ([5])]. 

For a mass profile with constant density ("tophat"), 
the force enhancement in Eq. (|52p is independent of ra- 
dius; for more general profiles however, this is not the 
case (see also [41[ for a detailed discussion). Fig. |4] shows 
the relative force enhancement ^(r) as function of ra- 
dius r/i?A in the case of an NFW halo, for the differ- 
ent DGP models ([H, HO], Tab. H]). We also show the 
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FIG. 4: ^DGp('') [Eq. (|52|l ] as function of the scaled radius 
r/i^A for an NFW halo, for the DGP models. The thin 
horizontal lines show the linearized expression ^DGP.iin = 
1 + 1/(3/3), while the thin lines deviating at small r show 
the results when using a "capped" density profile with reap = 
0.125 Mpc//i (see ^ IIIIBp . In all cases, we assumed A = 200 
and a concentration of ca = 5. For nDGP-1, we also show 
the effect of varying the concentration by ±20% (dotted lines). 
Note that ^dgp is indepedent of the halo mass. 



(r-independent) linearized value ^DGP.iin for the nDGP 
models which is recovered only at very large scales when 
the average density becomes < (for sDGP, ^oGP.iin ~ 
0.76 is beyond the range of the plot). 

Since ^DGP only depends on the average overdensity 
(5p/p, which is completely determined by A and ca, the 
force enhancement does not directly depend on the halo 
mass. Also, it is insensitive to the large-scale environ- 
ment of the halo. These arc two crucial distinctions from 
the /(i?) case. 

Fig. [5] shows J^.^. defined in Eq. (|27p as a func- 
tion of the overdensity A (keeping fixed at a value 
expected for a lO^^Mo/Zi halo), and the halo concen- 
tration. Clearly, J^.^ does depend somewhat on the 
halo profile and the overdensity criterion chosen. The 
general trend is that more concentrated halos lead to a 
stronger suppression of the modified forces, since they 
have higher average density at small radii. The same 
holds when increasing A. The dependence on c and A is 
strongest for nDGP-1 which also shows the strongest evo- 
lution of ^DGp('')- The dependence on the density profile 
has to be taken into account when comparing with simu- 
lation results fij llllBp . as well as for the comparison with 
observations which measure the dynamical mass within 
different (§|TV|. 
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FIG. 5: Averaged force deviation J^.^ [Eq. ((27|)] for DGP mod- 
els as function of the overdensity A {top panel) and the halo 
concentration C200 for an NFW halo {bottom panel). For the 
top panel we have scaled the concentration with Ra to keep 
Vs fixed (corresponding to C200 ~ 5.8). The thin hnes in the 
bottom panel show the results for a density profile capped at 
fcap = 0.125 Mpc//i for comparison with simulation results 
(assuming M200 = lO" M^/h; see ^ IIIIB|) . 



III. COMPARISON WITH SIMULATIONS 

In order to benchmark our theoretical expectations, we 
will now compare them to the results of the self-consistent 
N-body simulations of f{R) gravity presented in [itI l39j 
and of DGP [H, [l^l- For each model, we have simulated 
several box sizes. The number of runs for each model and 
box size, as well as the grid resolution are summarized in 
Tab. [Til Halos are identified using a spherical overden- 



TABLE 11: Number of runs for each box size and minimum 



mass cuts for cr,, and ' 



measurement. 





Model 


Lbox {h-^ Mpc) 


400 256 128 64 


#of 
runs 


sDGP 
nDGP' 


6 6 6 6 
6 6 6 6 
3 3 3 6 


Mh,min (10" Me/h) 


63.5 16.7 2.08 0.26 


rceii {h-^ Mpc) 


0.78 0.50 0.25 0.125 



"For each value of \fRo\. 

''For nDGP-1 and nDGP-2 each. 



sity halo finder as described in [32, |411 • The halo finder 
returns the center-of-mass position as well as mass Ma 
of the halo as determined from the particles within Ra, 
such that M I { Ait I^R\ ) = 'pls.. Our choice of A is the 
one adopted in [33, l4ll |: A = 300 for the f{R) simula- 
tions, and A = 200 in the DGP case. Our particle-mesh 
simulations are of limited resolution, and we can only use 
the best-resolved halos for our study, i.e. massive halos in 
the two smallest boxes. This limits our statistical sample 
of halos. 

First, in order to measure fj{r) as function of radius, 
we select the most well-resolved halos whose radii i? a are 
at least 10 grid cells, which is only satisfied for halos in 
our smallest box, Lbox = 64Mpc/ft,. For this box, this 
corresponds to a minimum mass of ^ 1.6 x 10"'^^Mq//i, 
which depending on the model results in a very small 
sample of 2-40 halos. For each halo, we then measure 
contributions to the potential energy W{r) in spherical 
shells around the centcr-of-mass via 



W{r) 



1 

TV 



^ (x» - x,0 • V^{yL^ 

-r|<Ar 



(54) 



where the sum runs over particles whose distance — 
|xi — x/i| from the ccnter-of-mass of the halo is within the 
radial bin, and N is the number of contributing particles. 
The derivative of the potential is evaluated at the posi- 
tion of each particle in the same way as it is done in the 
particle propagation of the N-body simulation (bilinear 
interpolation). In addition to W{r) derived from the dy- 
namical potential VP, we also measure the corresponding 
Newtonian quantity M^Ar(r), where the Newtonian po- 
tential VlfAT is determined from Eq. using the same 
density field. The ratio of the two is our estimated force 
enhancement: 



(r) 



W{r) 
WN{ry 



(55) 



To some extent, resolution effects can be expected to can- 
cel out in Eq. (j55p . Below we will show profiles down to 
f ~ fccW, though one should keep in mind that the ^ 
profiles cannot be considered reliable below r ^ 4rccii. 

Due to the resolution requirements and small sample 
size, we cannot study any evolution with mass in the 
^(r) profiles. The stringent resolution requirements can 
be relaxed somewhat if we only measure an average force 
enhancement, for example J^jj.- Assuming a scaling fol- 
lowing the virial theorem, we can either measure an aver- 
age of X- V^*, related to the trace of the potential energy 
tensor W (Eq. ([7])); or we can measure the velocity dis- 
persion, related to the kinetic energy given by Eq. ([9]). 
The first approach has the advantage that we can mea- 
sure J^jj. on a halo-by-halo basis, by calculating W using 
both the modified potential ^ and the Newtonian poten- 
tial ^jv (similarly to what was done for ^mcas (''))■ The 
estimator of J . . for a given halo is then defined by: 



' vir,meas 



(56) 
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FIG. 6: Velocity dispersion scaled to the virial expectation 
(ct^ oc M^^'^), measured for halos in the GR simulations, as a 
function of the halo radius in grid cells. Velocity dispersions 
are only reliably measured for the most well-resolved halos 
with i?3oo > 5.4 grid cells (indicated by the vertical line). 



where the sum runs over particles within the halo radius 
i?A- Note that the sum over particles automatically re- 
sults in a density weighting of g. Again, we expect that 
in this ratio resolution issues cancel to a certain extent. 
Some effects of the finite resolution will become apparent 
when comparing with the theoretical predictions below. 

The second approach, measuring the halo velocity dis- 
persions is also interesting since it gives an independent 
estimate of the modified forces. Specifically, we define 
the (one-dimensional) velocity dispersion of particles in 
a halo as follows: 



X-Xh|<flA 



"'^ = ^ 



(57) 
(58) 



X-X,j|<i?,^ 



where the sum runs over particles within the halo radius 
i?A, -/Vp is the number of those particles and V; — v/j 
denotes the velocity of the particle with respect to the 
center-of-mass of the halo. Note that in our normal- 
ization of (Tu, the kinetic energy Eq. (jH]) is given by 
T = 3/2 Missal- From the results of § Iff A[ we expect 
that when averaged over many halos, 



'i.,MG 



(59) 



where jyjQ is the velocity dispersion measured in the 



FIG. 7: gf(R){r) measured using Eqs. (|54[) - (|55p (thin lines), 
for the most well-resolved halos (-R300 > 10 grid cells) in the 
f{R) simulations, for the strong field |/ho| = 10~* and weak 
field I/ho I = 10^^. The mass range of the halos shown here is 
M300 = 1.6 - 7 X 10" Mq/H. The thick lines show the results 
of the relaxation code (for M300 = 3 x 10^"* Mq/H). 



modified gravity simulations, while qj^ is measured 
in the corresponding GR simulations. Note that in this 
measurement, we can only compare the average of many 
halos in the modified gravity simulations to that in GR, 
rather than calculating ^ on a halo-by-halo basis. Hence, 
Eq. (j59p results in a noisier measurement of J^j^, than 
Eq. (|56p . However, the particle velocity dispersion, which 
has gone through the relaxation and virialization process, 
is much more closely related to observables than the av- 
eraged gravitational force strength Eq. (j56p . which can 
never be measured directly in reality. Thus, it is worth- 
while to cross-check our results obtained from Eq. (|56p 
with the halo velocity dispersions. 

In order to determine for which halos we can reli- 
ably measure Eq. (f56|) and Eq. ([57| . we calculate the 
velocity dispersion of halos in the standard GR simu- 
lations and compare it to the expected virial scaling. 
Fig. ini shows the measured velocity dispersion scaled as 
a'^/AI'^^, as a function of the halo radius i?3oo in grid 
cells, for the different simulation boxes. Since the virial 
theorem has been found to hold in simulated halos [3ll |. 

a^/M^^ should be independent of the halo mass (S [H]). 
We see that, within the significant scatter, this is approx- 
imately true for halos that are sufficiently resolved. We 
thus place a radius cut of Ra > 5.4 grid cells. Since 
Ma = (47r/3) Api?^, this corresponds to a fixed mass 
cut for a given simulation box size, which is listed in 
Tab. m Fortunately, the statistics are then sufficient to 
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FIG. 8: Same as Fig. [T] but for the intermediate field value 
|/iio| = 10"^. We iiave separated the halo sample into lower 
mass halos with M300 ~ 1.6 — 2. 5x 10^"* Mq/H, and two higher 
mass halos with M300 = 6 - 7 x 10^'' Mq /H. 

study J^.^. and as functions of mass. 

A. f(R) 

We begin with the measurement of ^/(fl) for the well- 
resolved halos. Fig.[7]shows the simulation measurements 
and predictions of the spherical relaxation code, for the 
strong field {\fm\ = 10"'*) and the weak field (10"^^). 
As expected from Fig. [3l the halos in this mass range 
(-^^300 ~ 1.6 — 7 X Mq/K) are unscreened in the 
strong field simulations, gf(R) = 4/3, and screened in the 
weak field, sf(R) — > 1. In the latter case, there is a regime 
around 1-3 i?3oo where the screening is not complete. At 
larger distances, the Yukawa suppression again becomes 
noticeable (Ac ~ 3 Mpc for \fm\ = 10~^). The numeri- 
cal results for the spherically symmetric case match the 
overall behavior well for both field values, although there 
is a hint that it slightly overestimates the screening in 
the weak field case. 

The results for the intermediate field value \fm\ ~ 
10~^ are shown in Fig. [51 this case is most interesting 
since the fewx 10"'^"' Mq/Zi halos are in the transition re- 
gion between the screened and unscreened limits (Fig.[3|). 
Hence we have split the halo sample into a lower mass 
sample around 2 x IO^'^Mq/Zi and a high mass sample 
with two halos around 7 x IQ^'^ Mq/H. Clearly, the scat- 
ter in the modified force profiles is significant. Neverthe- 
less, the stronger screening effect in the higher mass halos 
is noticeable. The spherical relaxation results (§ III Cp . 




10'" 1015 
M300 [Mo/h] 



FIG. 9: J^jj, measured via Eq. (|56|l for well-resolved halos 
(-R300 > 5.4 grid cells) in the ,f{R) simulations (points). The 
results confirm the theoretical predictions from § III CI shown 
as lines (spherical relaxation results from Fig. [S)). The cir- 
cled points are halos which have a more massive halo in their 
immediate vicinity (see text). 

which were calculated separately for the mean halo mass 
of each sample, match the full simulation profiles remark- 
ably well. At small radii, the transition to the fully 
screened values is apparently too steep. A possible expla- 
nation for this is that the halos in the N-body simulations 
are not truly spherical, but in general triaxial. A triax- 
ial halo will have a somewhat shallower potential well, 
reducing the chameleon screening effect. Furthermore, 
the screening will happen at different radii along the dif- 
ferent axes, so that a potentially sharp transition in the 
spherical case is washed out over a certain radius range. 
In addition, the innermost Newtonian potential well is 
not as deep in the simulations as predicted for a perfect 
NFW halo due to the finite resolution. We also reiterate 
that the profiles only become reliable at r ~ 0.3 — 0.4i?3oo 
for these halos. 

Next, we look at 7 . in the larger sample of halos. 

' '-'vir,mcas ^ ^ 

Fig. [9] shows the results for the three field values, and the 
predictions oij^.^^ ^^.^^(M) from the spherical relaxation 
code. We again see a very good match for all field val- 
ues and over the entire mass range probed by this halo 
sample, 3 x Mq/H < M300 < 3 x lO^'^Mo/ft,. For 
the intermediate field value, which again shows the most 
interesting behavior in this mass range, we see that the 
screening effect is slightly overpredicted in the spherically 
symmetric approximation. Again, this could be due to 
halo triaxiality and to resolution effects which reduce the 
value of 4* AT in the inner parts of the halo. 
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FIG. 10: Scaled velocity dispersion cr^/M^''^ measured in GR 
and f{R) simulations. The measurements were scaled to val- 
ues expected for NFW halos (§ HT]): erg = 1.79 x 10"^c^ 
Mo = 10^^ Mq/H. The dotted black line shows a constant 
fit to the GR results. Solid and dashed lines show the predic- 
tions of the model of III CI scaled by the GR value. 



For the intermediate field, some outliers are seen in 
Fig. ini These halos, especially around 3 — 8 x 10^^ Mq/H 
show a stronger screening of the modified forces than the 
large majority of halos at that mass, and stronger than 
predicted for isolated spherical NFW halos. This would 
seem consistent with halos being screened by a larger 
scale potential well in which they are situated. To test 
this hypothesis, we have selected halos in the intermedi- 
ate field simulations which have a more massive neigh- 
boring halo in their immediate vicinity. More precisely, 
we ask that 



i?A. 



N 



RA,h 



< 1, 



(60) 



where h denotes the halo itself and N denotes the closest 
neighboring halo^ with Ma.n > Maji- These halos make 
up less than 5% of the whole sample and are circled in 
Fig. [9l In fact, three of the outliers have a close massive 
neighbor, which is strong evidence for the hypothesis of 
environmental effects as cause for the enhanced screening 
(the fourth most obvious outlier has djy ~ 1-2)- 



FIG. 11: Same as Fig. [TOl but for the intermediate f(R) field 
value |/flo| = 10-^ 



Finally, we can look at the effect of the modified forces 
on the particle velocity dispersion of halos, a noisier mea- 
surement but one that probes the effect after the repro- 
cessing through gravitational collapse and virialization. 
Fig. [TUl shows the scaled velocity dispersion a^/Ad^^^, for 
the same halos as in Fig. |9] and scaled to values expected 
from § ini as a function of mass. We show the results 
for GR simulations as well as the weak (|/iio| = 10^^) 
and strong field (|/i?o| = 10^'') f{R) cases. After fit- 
ting a constant to the GR simulations, we multiply the 
theoretical predictions by this constant. Albeit noisy, 
the results of the halo-by-halo measurement of J^.^ ^^^j 
are confirmed: for the strong field, all halos are in the 
linearized field regime where forces, and hence tr^ are en- 
hanced by a factor of 4/3. For the weak field, all halos 
except at the very lowest masses probed by the simu- 
lations are in the chameleon regime. In case of the in- 
termediate field (|/i?o| = 10~^), Fig. [TT] shows that the 
transition between screened and unscreened regimes at 
fcwx 10^^ A/q/Zi is indeed seen in the halo velocity dis- 
persions as well. These results confirm that the theoret- 
ical predictions for the modified gravitational force can 
in principle be probed by observable quantities such as 
velocity dispersions (see ? IIV[) . 



^ This criterion formally says that the halos are overlapping. Such 
an overlap is unavoidable wrhen defining halos via spherical over- 
densities. In our halo finding algorithm, the particles in the 
overlap region are not double-counted, but counted towards the 
more massive halo. 



B. DGP 

The force modifications in DGP are, to first order, in- 
dependent of the halo mass and environment. However, 
they do depend on the detailed halo profile. Before com- 
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FIG. 12: ^DGp{r) measured using Eqs. ((54|) - (f55l) for the 
most well-resolved halos (i?2oo > 10 grid cells) in the DGP 
simulations (thin lines). The thick lines show the predic- 
tion of Eq. using a capped NFW profile with reap = 
0.125 Mpc//i (see text). The thin horizontal lines show 
.jDGP.iin for each model. The halos shown here have masses 
M2oo'= 1.6 - 7 X 10" Me/h. 



paring the predictions from ? Ill Dl with the simulation re- 
sults, we have to take into account the effects of the finite 
resolution. While the NFW profile we used throughout 
§|TT]is a very good match to high-resolution simulations, 
in our fixed-grid simulations the density profile is in fact 
softened on scales of a grid cell. This softening of the den- 
sity profile will affect ^dgp through the average over- 
density within r. Thus, for comparison with the simula- 
tion results we assume a "capped" density profile instead 
of NFW all the way to r = 0. More precisely, we cap the 
density profile at a constant value of pcap = PNFw(^cap) 
for r < Tcap^- For the halos measured in the smallest box, 
a natural choice is reap = ''ccU = 0.125 Mpc/ft, (Tab. HB. 
Fig. |4] shows the effect of this softened density profile on 
^DGp('')- In particular, it increases the force modifica- 
tion since the inner density is suppressed, thus artificially 
weakening the Vainshtein mechanism. 

Fig. [T^] shows the measured ^dgp (f) from the simu- 
lations, together with the predictions using the capped 
density profile. First, it is evident that the scatter in ^ is 
much smaller in DGP than it is for f{R), due to the lo- 



^ Note that the halo radius for a given mass is slightly increased 
when using the capped density profile, in order to match the 
fixed overdensity A = 200. 
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1 _ '■ nDGP-2 
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FIG. 13: ^^j^ measured using Eq. (|56|) for well-resolved 
halos (-R200 > 5.4 grid cells) in the DGP simulations (points). 
The shaded bands show the model predictions from ij lll Dl with 
a variation in the halo concentration by ±20%. We assumed 
capped NFW profiles with reap = 0.125 Mpc/ft. 

cality of the Vainshtein mechanism. For all three models, 
the agreement of the simple spherically symmetric NFW 
model with the simulations is impressive. Note that we 
have not adjusted any parameters to match the simu- 
lation results; this measurement thus also constitutes a 
nontrivial test of the DGP simulations. At r '--^ -R200! the 
theoretical prediction slightly underestimates the sup- 
pression of the force modification (by 1 — 3%), which is 
presumably due to slight differences in the actual density 
profiles from the one assumed in the predictions (pure 
spherical NFW profile). The large scatter at r > 2i?2oo 
is due to the effect of gravitationally unbound ambient 
matter in the environment of the halos, which dominates 
Sp at these distances. Note that in particular for uDGP- 
1, the Vainshtein mechanism does not completely sup- 
press the force modifications within halos even on scales 
as small as ^ 100 kpc. 

We now turn to J^^^. as measured from Eq. (j56p in the 
halo sample with i?2oo > 5.4 grid cells. Fig. [13] shows 
the measurements for the three DGP models. As ex- 
pected, J^.^ j-|Qp is approximately constant as a function 
of mass. The model predictions from § IIIDI are shown 
as gray bands. Here, we have used the concentration re- 
lation Eq. p6|) (more precisely, c = max{4, c(A'/)}), and 
the width of the band reflects a ±20% spread in con- 
centration. We again assumed a capped NFW profile 
with reap = 0.125 Mpc//i, the effects of which are no- 
ticeable as a slight increasing trend of J^.^(M) in going 
towards the low- mass end for uDGP-l. Note that here 
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FIG. 14: Scaled velocity dispersion al/M'^^'^ measured in 
GR and nDGP-1 simulations (cro and Mo are as defined in 
Fig. nop. The black dotted line shows the average value mea- 
sured for the GR simulations. The red line shows this value 
multiplied by J^.^ oGPiin = 1 + 1/(3/3). The shaded band 
shows the corresponding prediction for J^.^ j^^p from Fig. [131 



we have included halos from different simulation box sizes 
Lhox — 64 — 256 Mpc//i. though the majority comes from 
the smallest box. Hence, one might wonder whether dif- 
ferent values of reap are required for different box sizes. 
However, within the hmited statistics the measurements 
of 7 . from halos in different box sizes are in agreement. 

^ vir ^ 

Hence, the data do not seem to require such a correc- 
tion. We conclude that, within the uncertainties due to 
the halo density profiles, the measured values of J^.^ 
are entirely consistent with the predictions. Furthermore, 
the scatter in the measured J^j^, appears consistent with 
that expected for intrinsic variations of halo density pro- 
files (Ac/c- 0.2). 

The results for J^.^. j^^^p arc confirmed by the particle 
velocity dispersions of halos. Fig. [14] shows the scaled 
velocity dispersion cr^/A/^/"^ in nDGP-1 and the corre- 
sponding GR simulations. For comparison, we also show 
the result for the linearized DGP simulations, which use 
the scale-independent (but redshift-dependent) force en- 
hancement obtained when linearizing the DGP equations 
[isl [20| . Within the significant scatter in the measure- 
ment, we found no significant evolution of the force en- 
hancement with mass, as expected given the small trends 
with mass in Fig. 1131 

In order to quantitatively compare the simulation re- 
sults with model predictions, we determined the mean 
of a^jM'^l'^ for each simulation type. In each case, the 



error on this mean is obtained by dividing the RMS scat- 
ter by V-^haios- The measured ratio of the scaled velocity 
dispersion in the DGP simulations to that in the GR sim- 
ulations is found to be 



. .(full DGP, meas) = 1.212 ± 0.014. 



(61) 



This is indeed close to the range of the theoretical predic- 
tions (for a capped NFW profile), 1.18-1.2 (Fig.[T31). As 
expected; the ratio measured in the linearized DGP sim- 
ulations, J^^j^, (lin. DGP, meas) = 1.288± 0.014 is in excel- 
lent agreement with the predicted value of 1 -I- 1/(3/3) = 
1.276. Similar conclusions hold for the velocity disper- 
sions measured in the sDGP and nDGP-2 simulations, al- 
though the results are less constraining due to the smaller 
force modifications |7 . — 1| in those models. 



IV. APPLICATION TO OBSERVATIONS 

Observables linked to dynamical masses can be broadly 
classified into two categories. First, one can measure 
the velocity distribution of collisionlcss "tracer particles" , 
such as galaxies within galaxy clusters, or stars within 
galaxies. For a dynamically relaxed system, the kinetic 
energy T inferred from the velocity distribution is pro- 
portional to the potential energy W (§ III A[) . which can 
be converted into a mass estimate MA.dyn (we again as- 
sume a mass definition in terms of an average interior 
density pA). Several assumptions have to be made in 
order to obtain the mass estimate. First, one has to as- 
sume the galaxies or stars are unbiased tracers of the 
full matter velocity field (including dark matter). Since 
member galaxies of a cluster generally reside in overdense 
substructure (subhalos) of the cluster halo, their veloci- 
ties might differ systematically from that of the overall 
matter. Simulation studies [la, 113] have shown that this 
velocity bias is expected to be on the order of ^ 10% 
or less, depending on how galaxies are selected. Further, 
one has to make assumptions about the density profile 
shape, and the anisotropy of the velocity distribution, 
since only the linc-of-sight component of the velocity is 
observed. Nevertheless, our idealized measurement of the 
dark matter velocity dispersion in the simulations shows 
that at least in principle, is indeed a good tracer of 
the modified force 1 . . 

o vir 

Another set of observations linked to dynamical masses 
is measurements of the hot ionized gas in galaxy clusters. 
One technique is to detect the thermal bremsstrahlung in 
X-rays; another is to measure the upscattering of CMB 
photons off the hot electrons via the Sunyaev-Zeldovich 
(SZ) effect. In both techniques, one measures a line- 
of-sight integral of the electron pressure, with an ad- 
ditional weighting by the electron density in the case 
of X-rays (since the rate of bremsstrahlung emission is 
(X Tie Up = rig). With some assumptions on the density 
profile for the baryons. X-ray and SZ signals can be con- 
verted into a measurement of the electron pressure as 
function of r. Instead of the virial theorem which holds 
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for a collisionless system, we now use hydrostatic equilib- 
rium which is a good assumption at least for dynamically 
relaxed systems: 



dP 

dr 



Pgas 



dr ' 



(62) 



where P is the total pressure and pgas is the mass den- 
sity of the gas, respectively. The difficulty observation- 
ally is in measuring the left-hand side of Eq. (|62|) : only 
the thermal contribution to P, Pthcrm ~ ng^skT is di- 
rectly measurable, while non-thermal contributions from 
e.g. cosmic rays, bulk flows, and magnetic fields are much 
harder to estimate. Nevertheless, with appropriate sys- 
tematic error bars, Eq. (p^ is a probe of the gravitational 
force /dr. 

In summary, a variety of observations lead to estimates 
of certain weighted averages of the gravitational force. 



obs 



j d3xpobs(x)x- V*(x), 



(63) 



where pobs is an effective weight function. In case of 
X-ray and SZ measurements, it is related to p^^^ and 
Pgas, respectively, but will be modified by instrumental 
effects such as the limited instrument aperture. Simi- 
larly, for galaxy velocity dispersions in clusters, pobs is 
proportional to the number density of observed galaxies 
(again, with observational weights and boundaries folded 
in). 

Now we can use Eq. pU)) together with the fact that 
cx M|^/^ so that W a M^'^ [Eq. ©]. Then, if the 
observational mass estimate is done based on ordinary 
gravity, so that in GR the mass estimate equals the true 
mass Ma, the resulting mass estimate MA,dyn in modi- 
fied gravity is in fact 



Ma 



dyn - ^^A- 



(64) 



Here J^^^ is a weighted integral using Eq. pT|) . with p 
replaced by the effective weight pobs- Note that in general 
^obs ^^^^ depend on the true mass Ma itself. 

Since the true mass can in principle be obtained from 
weak or strong lensing, a comparison of lensing mass with 
the dynamical mass Eq. (|64| can be used to measure 
the modified forces in f{R) or DGP. Again, it is impor- 
tant to take into account the unavoidable observational 
weighting that is being done in the measurements of both 
MA,dyn and Ma- 

Recently, the SLAGS sample of elliptical galaxies act- 
ing as strong lens es 1481 h as been used to constrain de- 
viations from GR [28l l49l| . Furthermore, using a similar 
argument as the thin shell condition Eq. ()36p G\ . [49l | 
has shown that these measurements constrain the f{R) 
model considered here at the level of \fBa\ ^ 2 x 10^^. 
For these constraints one has to make some assumptions 
on the potential well of the lens galaxy, for example that 
it is dominated by the density distribution of the inner 
few kpc, thus neglecting any larger scale potential well. 



As we have seen, the magnitude of the force modifica- 
tion in /(i?) can depend somewhat on the environment. 
In particular, we found that a subset of halos around 
3— 8 X 10^"^ Mq jh (at the low mass end of the range acces- 
sible to the simulations) is screened much more strongly 
than expected for isolated halos, consistent with an effect 
of the large scale environment. Nevertheless, strong lens 
galaxies offer a quite powerful probe of gravity on kpc 
scales, if the environmental effects can be understood. 

On larger scales, the comparison of dynamical and 
lensing masses of massive galaxy clusters can be inter- 
esting since they dominate their local environment, so 
that environmental effects should be negligible. Also, for 
cluster-scale masses we were able to validate our theo- 
retical models for J^.^ directly with the modified gravity 
simulations (SHU). However, for clusters it is preferable 
to measure the dynamics and lensing at large scales: first, 
the deviations from GR quickly shrink close to the clus- 
ter core owing to the chameleon and Vainshtein mech- 
anisms; second, baryonic effects on the observables and 
the density profile, such as cooling and AGN feedback, 
are expected to be less significant at greater distances 
from the cluster center. 

It is also possible to use dynamic mass estimates of 
clusters by themselves, without direct comparison to 
lensing masses. As shown in j50l - l55| . the abundance of 
massive clusters is a sensitive probe of the growth of 
structure as well as gravity. When comparing the ob- 
served cluster mass function measured using a dynami- 
cal mass measure with modified gravity predictions, it is 
necessary to take into account the effect of the modified 
forces on the mass estimates as well. In order to estimate 
the effect on the observed mass function, we use Eq. (|64p , 



setting , 



the idealized quantity we have mod- 



eled and calibrated with simulations. Dynamical mass 
measures (i.e. velocity dispersions) in our simulations 
are noisy (§ IIII[) . thus we have simply rescaled the mass 
of each halo in the modified gravity simulations by our 
theoretical model of J . (M) for the given cosmology. 

Fig. [15] shows the relative enhancement of the mass 
function in /(i?) gravity with respect to ACDM, when 
measured using lensing masses (i.e. "true" M300) and dy- 
namical masses, for'' |/ijol = 10~^ and 10~^. Glearly, the 
observed abundance of halos is further enhanced when 
measured in terms of dynamical masses. In the mass 
range where halos are unscreened, the mass function en- 
hancement is boosted by a factor of two or more. This 
is because the dynamical mass estimate is a factor of 
(4/3)^/^ « 1.19 higher than the lensing mass in f{R) 
gravity in the unscreened case, in conjunction with the 
steeply falling mass function. Gonstraints on f{R) grav- 
ity from X-ray clusters could thus be significantly im- 



^ Since all halos above 10^'' Mq /h are screened for the small 
field |/i{o| = 10~^, the dynamical mass function is essentially 
equal to the lensing mass function for most of the mass range, 
and is not repeated here (see |32|| ). 
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FIG. 15: Mass function enhancement in f{R) relative 
to a ACDM cosmology with the same expansion history, 
ninA/(/(i?)) /ninA/(ACDM) - 1, for \fm\ = 10"* {top panel) 
and 10~^ {bottom panel). The points show simulation results, 
while the shaded bands show spherical collapse predictions 
[3^ (see text). Results are shown for the mass function nin 
in terms of the lensing mass, and nin ^yn terms of the 
dynamical mass (see text; A = 300 in both cases). 

proved by using the dynamical mass function instead of 
the true or lensing mass function which was used in js^ ] . 
Note the sharp turnover in the mass function enhance- 
ment for the intermediate field value. This transition 
due to the onset of the chameleon mechanism is already 
present in the lensing mass function [s^ . Since J^^^ tran- 
sitions from 4/3 to 1 in this mass range as well, the effect 
is enhanced in the dynamical mass function. 

The shaded bands in Fig. [T^ show the spherical col- 
lapse predictions presented in |32| . These are based on 
the linear f{R) matter power spectrum together with the 
Sheth-Tormen prescription, using two sets of collapse pa- 
rameters derived for limiting cases of spherical collapse 
in f{R) (enhanced forces throughout, and unmodified 
forces). We rescaled the predictions in terms of lensing 
mass given in [3^ to the dynamical mass via 

_ dn din Ma 

As expected, the predictions in terms of dynamical mass 
perform equally well as those for the lensing mass. Since 
our prediction for J^.^ includes the chameleon mecha- 
nism, the predictions for the intermediate field value 
show a corresponding transition at approximately the 
right mass. Still, the predictions do not match the simu- 
lation results completely due to the shortcomings of our 
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FIG. 16: Same as Fig. [TS] but for the DGP models nDGP- 
1 {top panel) and nDGP-2 {bottom panel). The simulation 
results for the dynamical mass (red triangles) have been dis- 
placed horizontally for clarity The shaded band shows the 
spherical collapse model of [4l| . 

simple spherical collapse model f32j . 

Fig. [in] shows the corresponding results for the two 
normal-branch DGP models nDGP-1 and nDGP-2. The 
effect is less dramatic on the DGP mass function, since 
J^jj. in DGP is generally smaller than in f{R). Never- 
theless, the impact especially for nDGP-1 is significant, 
an additional implying an abundance boost of ^--^50% at 
high masses. The shaded bands in Fig. [12] again show a 
spherical collapse model [4l|, which uses the analytical 
solution for the modified forces in DGP in the spherically 
symmetric case as one limiting case of spherical collapse 
in DGP. The other limit is given by using the linearized 
expression for the modified forces. Again, the spherical 
collapse model performs equally well for the mass func- 
tion in terms of dynamical mass as for the lensing mass 
function. 



V. CONCLUSIONS 

In this paper, we have studied the dynamics of matter 
within bound cosmic structures, i.e. dark matter ha- 
los, in f{R) and DGP. The potential governing mat- 
ter dynamics can differ from the lensing potential by 
20 — 30% in these models. These unique signatures of 
modified gravity can be observed by comparing dynam- 
ical and lensing mass estimates of clusters or galaxies. 
Furthermore, they strongly influence the observed abun- 
dance of massive clusters when measured via dynamical 
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mass proxies such as X-rays or the SZ effect. For exam- 
ple, the enhancement of the cluster abundance in f{R) 
(with respect to ACDM) at a fixed dynamical mass can 
be roughly twice that measured if the mass is based on 
lensing measurements. These signatures in the dynamics 
are also relevant for large-scale structure observations, 
such as the redshift-spacc power spectrum or correlation 
function on small scales. 

However, since halos are highly non-linear objects, the 
peculiar chameleon and Vainshtein mechanisms play a 
crucial role, as they are necessary in order to restore Gen- 
eral Relativity in high-density environments. Thus, the 
dynamics in these models can only be rigorously studied 
through N-body simulations which include the non-linear 
mechanisms of f{R) and DGP consistently. 

In the case of f{R), the chameleon mechanism is trig- 
gered once the depth of the potential well is comparable 
to the background value of the scalar field. The suppres- 
sion of the force modifications within a halo thus depends 
not only on the halo mass but also its environment. Con- 
sequently, we found significant scatter from halo to halo 
in the force modification ^ measured in the f{R) simula- 
tions. Furthermore, we identified a subset of halos which 
are in the close vicinity of massive neighbors, and which 
show a much stronger suppression of the force modifica- 
tions than expected for isolated halos. In the majority 
of cases however, the simulation results confirm the basic 
expectation that halos are "unscreened" below a certain 
threshold mass determined by the potential well and the 
field value, whereas GR is restored at higher masses. Fur- 
thermore, a simple model based on the spherically sym- 
metric solution of the field equations provides a good 
match to the scale- as well as mass-dependence of the 
force modifications in f{R). 

In DGP, the non-linear suppression of the force modi- 
fications through the Vainshtein mechanism is much less 
dependent on halo mass and details of the large scale en- 
vironment. Instead, the crucial quantity is the average 
mass density within a given radius. Thus, uncertainties 
in the semi-analytic predictions for DGP are mainly due 



to the density profile, and are already quite small. When 
taking into account the force resolution of the simula- 
tions, our predictions provide an excellent fit to the sim- 
ulation measurements. Since the basic assumptions of 
the model, in particular spherical symmetry seem to hold 
well, we expect that force modifications can be predicted 
very accurately in DGP, provided the density profile is 
known sufficiently well. 

Given that our semi-analytic models appear to capture 
the mass- and scale-dependence of the modified forces 
correctly for both f{R) and DGP, they can be useful in 
extending predictions beyond the limits of resolution and 
parameter space of the simulations. This will be neces- 
sary in particular for the comparison with observations. 

While this study is specific to f{R) and DGP, it shows 
the qualitative features expected in observations of dy- 
namics from viable modified gravity models, which em- 
ploy a non-linear mechanism to restore GR locally. In 
the outer regions of massive clusters, as well as in lower 
mass objects, these models generally predict order unity 
deviations from GR. Observations in this regime thus of- 
fer the perspective of closing the last remaining loopholes 
for significant modifications to gravity on large scales. 
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